#include<stdio.h>
#include<stdlib.h>
#include <mpi.h>
#include <omp.h>
#include <string.h>

#define BLOCK_LOW(id,p,n) (((id)*(n))/(p))
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,(p),(n)) - 1)
#define NTHREAD 2

float   **Alocar_matriz_real (int, int);
float   **Liberar_matriz_real (int, int, float **);
int	LeituraArquivo(int,int,int,FILE *,int,float **,float **,int,int);
float	*Alocar_vetor_real(int);
int	*Alocar_vetor_int(int);
float	*Liberar_vetor_real(int, float *);
int	*Liberar_vetor_int(int, int *);
float	Numero_maximo(int ,float *);

int main(int argc, char *argv[])
{

        int     ordemMatriz;                    // J_ORDER
        int     filaAvaliacao;                  // J_ROW_TEST
        float   erro,mr;                        // J_ERROR mr= erro por inreracao
        int     interacoes;                     // J_ITE_MAX
        int     nl;                             // Indica o numero de  linha que cada processo deve ler do arquivo
        float   **matriz,**matOriginal;         //matriz que contem os dados
        int     i,j,k,posicao;                  //manipulacao loops
        float   aux;
        int     daux;
	int	process_fila;			//processo q contem a fila de avaliacao
        FILE    *fp;
        float   *x;                             //resultado local das iterações
        float   *xg,*xgAntes,*dif;              //resultado global das iterações e criterio de parada
        int   flag = 0;
        int     *rcv_disp;                      //indices do MPI_Allgatherv
        int     *rcv_cnt;  
	int	resto;     

        /*Variáveis MPI */   
        int             dest,noProcesses, processId, src, tag;
        MPI_Status      status;   

        /* Inicializacção do MPI */   
        MPI_Init(&argc, &argv);   
        MPI_Comm_size(MPI_COMM_WORLD, &noProcesses);
        MPI_Comm_rank(MPI_COMM_WORLD, &processId);

        //abre o arquivo para leitura - todos os nos possuem esse arquivo
        if ( (fp = fopen(argv[1],"r")) == NULL ) 
        {
                printf("Erro ao abrir o arquivo.\n");
                exit(1);
        }

        //leitura dos dados iniciais - Todos os nos precisao desse valor
        fscanf(fp,"%d",&ordemMatriz);
        fscanf(fp,"%d",&filaAvaliacao);
        fscanf(fp,"%f",&erro);
        fscanf(fp,"%d",&interacoes);
        
	resto = ordemMatriz%noProcesses;

        //divide o numero de linha que cada processo deve ler
        nl=(int)(ordemMatriz/noProcesses);
        if(resto>0) {
                        flag = resto;
                        if(processId==0) {
                                flag = 0;
                                nl+=resto;
                        }
                                
        }
	
	// define a posicao inicial de cada nó no vetor de resultados, onde cada nó começa
	posicao=(int)BLOCK_LOW(processId,noProcesses,ordemMatriz-flag)+flag;
                
        //reserva memoria necessaria para cada nó 
        matriz=Alocar_matriz_real (nl + 1,ordemMatriz); // +1 porque tb contera os dados da matriz
	matOriginal=Alocar_matriz_real (nl + 1,ordemMatriz); // temos essa matriz apenas para manter 
							     // os dados inicias, assim dando a saida desejada
        x = Alocar_vetor_real(nl);
        xg = Alocar_vetor_real(ordemMatriz);
	xgAntes = Alocar_vetor_real(ordemMatriz);
	dif= Alocar_vetor_real(ordemMatriz);
        rcv_disp = Alocar_vetor_int(noProcesses);
        rcv_cnt = Alocar_vetor_int(noProcesses);
        
        //define os vetores utilizados em MPI_Allgetherv
        if (resto>0) {
                rcv_disp[0] = 0;
                if (processId == 0) rcv_cnt[0] = nl;
                else rcv_cnt[0] = nl+resto;
                for ( i = 1; i < noProcesses; i++ ) {

                        rcv_disp[i] = (int)BLOCK_LOW(i,noProcesses,ordemMatriz-resto)+resto;
                        if (processId == 0) rcv_cnt[i] = nl-resto;
                        else rcv_cnt[i] = nl;

                }
        }
        else {

                for ( i = 0; i < noProcesses; i++ ) {

                        rcv_disp[i] = (int)BLOCK_LOW(i,noProcesses,ordemMatriz);
                        rcv_cnt[i] = nl;

                }

        }


        //Leitura arquivo
        process_fila=LeituraArquivo(processId,noProcesses,nl,fp,ordemMatriz,matriz,matOriginal,filaAvaliacao,flag);

        //divide tudo pela diagonal principal
        for(i=0;i<nl;i++)
        {
                aux=matriz[i][posicao + i];
                for(j=0;j<ordemMatriz;j++)
                        matriz[i][j]=matriz[i][j]/aux;

                matriz[nl][posicao + i]/=aux;
        }
        

        //na primeira interacao, apenas atribuimos a x o valor do vetor
	k = 1;//primeira interacao
        for(i=0;i<nl;i++) 
                x[i] = matriz[nl][(int)BLOCK_LOW(processId,noProcesses,ordemMatriz-flag)+i+flag];


        //enche o buffer para a proxima interaco
        MPI_Allgatherv(x,nl,MPI_FLOAT,xg,rcv_cnt,rcv_disp, MPI_FLOAT,MPI_COMM_WORLD);

	omp_set_num_threads(NTHREAD);
	
	while(k<interacoes)
	{
	 
			#pragma omp parallel for shared(xg,ordemMatriz) private(i) firstprivate(xgAntes)
			for(i=0;i<ordemMatriz;i++)
				xgAntes[i]=xg[i];


			#pragma omp parallel for shared(ordemMatriz, posicao,matriz,nl) private (i,j) firstprivate(aux,x)
			for(i=0;i<nl;i++)
			{
				aux=0;
				//aqui cabe o paralelo for
				for(j=0;j<ordemMatriz;j++)
					if((posicao + i) != j)
						aux=aux+matriz[i][j]*xg[j];

				//Xdepois(i)=Vetor(i) - Soma(j=0 -> n)de{matriz[i,j]*Xantes(j)} para i<>j
				x[i] = matriz[nl][posicao + i]-aux;	
			}

		k++;
		//carrega o buffer
		MPI_Allgatherv(x,nl,MPI_FLOAT,xg,rcv_cnt,rcv_disp, MPI_FLOAT,MPI_COMM_WORLD);
		
		//criterio de parada
		
		#pragma omp parallel for shared(ordemMatriz,xg,xgAntes) private(i) firstprivate(dif)
		for(i=0;i<ordemMatriz;i++)
		{
			dif[i]=xg[i]-xgAntes[i];
			if (dif[i] <= 0) dif[i] = (dif[i] * -1); //abs(dif[i])
		}


		
		mr=Numero_maximo(ordemMatriz,dif)/Numero_maximo(ordemMatriz,xg);
		if (mr <= 0) mr = (mr * -1); //abs(mr)
		if(mr<=erro)
			break;
		//printf("Iterations: %f -- ",mr);
	

	}
	
       	//impressao resultados
	if(processId==process_fila)//apenas um nó imprime os resultados
	{
		printf("Iterations: %d\n",k);
		aux=0; 
		if(processId==0)
			j=filaAvaliacao-(nl*process_fila)-1;
		else
			j=filaAvaliacao-(nl*process_fila)-1+resto;

		for(i=0;i<ordemMatriz;i++)
			aux=aux+matOriginal[j][i]*xg[i];

		printf("Iterations: %d\n",k);
		printf("RowTest: %d => [%f] =? %f\n\n",filaAvaliacao,aux,matOriginal[nl][filaAvaliacao-1]);
	}

        matriz=Liberar_matriz_real (nl+1,ordemMatriz,matriz);
        x=Liberar_vetor_real(ordemMatriz,x);

        MPI_Finalize( );
	fflush(stdout);
        
        return 0;
}

void atualizar_vetor(float *xg, float *x, int pos, int nl, float **matriz) {

        int i = 0;

        for ( i = 0; i < nl;  i++ ) {           

                x[i] = matriz[nl][pos];

        }

}

int LeituraArquivo(int processId,int noProcesses,int nl,FILE *fp,int ordem,float **mat,float **m,int filaAvaliacao, int flag)
{
        int tamanho;
        int i,j,local;
        long int aux;
        local=-1;
        if(processId == 0)
	{
                for(i=0;i<nl;i++)
                        for(j=0;j<ordem;j++)
                        {
                                fscanf(fp,"%f",&mat[i][j]);
				m[i][j]=mat[i][j];
                        }
		if(filaAvaliacao <= i)
		{
			local=0;

		}
	}
        else
        {

                tamanho=2*ordem-1; //no minimo, cada linha tera um digito segido de um espaco, o ultimo digito nao tera espaco (-1)
                //posiciona na posicao correta          
                for(i=0;i<((processId*nl)+flag);i++)
                {
                        fseek(fp, ftell(fp)+tamanho, SEEK_SET);
                        while (getc(fp)!= '\n');
                }
		
                //agora basta ler o arquivo
                for(i=0;i<nl;i++)
                        for(j=0;j<ordem;j++)
                        {
                                fscanf(fp,"%f",&mat[i][j]);
				m[i][j]=mat[i][j];
                        }
		if((filaAvaliacao > ((processId*nl)+flag)) && (filaAvaliacao <= ((processId*nl)+flag+i)))
		{
			local=processId;
		}
        }
        //vai no final do arquivo e le de traz pra frente o vetor final
        fseek(fp, 0, SEEK_END);
        fseek(fp, ftell(fp)-2, SEEK_SET);       
        
        for(j=ordem-1;j>=0;j--)
        {
                while (getc(fp)!= '\n')
                        fseek(fp,ftell(fp)-2, SEEK_SET);

                aux=ftell(fp);
                fscanf(fp,"%f",&mat[nl][j]);
		m[nl][j]=mat[nl][j];
                fseek(fp, aux-2, SEEK_SET);     
        }

	return local;
}

float **Alocar_matriz_real (int m, int n)
{
        float   **v;  // ponteiro para a matriz 
        int     i;    // variavel auxiliar    

        if (m < 1 || n < 1) // verifica parametros recebidos
        { 
             printf ("** Erro: Parametro invalido **\n");
             return (NULL);
        }

        // aloca as linhas da matriz 
        v = (float **) calloc (m, sizeof(float *));
        if (v == NULL) 
        {
             printf ("** Erro: Memoria Insuficiente **");
             return (NULL);
        }
        
        // aloca as colunas da matriz
        for ( i = 0; i < m; i++ ) 
        {
                v[i] = (float*) calloc (n, sizeof(float));
                if (v[i] == NULL) 
                {
                        printf ("** Erro: Memoria Insuficiente **");
                        return (NULL);
                }
        }
        return (v); // retorna o ponteiro para a matriz
}

float **Liberar_matriz_real (int m, int n, float **v)
{
        int     i;  // variavel auxiliar

        if (v == NULL) return (NULL);
        if (m < 1 || n < 1) // verifica parametros recebidos
        {  
                printf ("** Erro: Parametro invalido **\n");
                return (v);
        }
        for (i=0; i<m; i++) free (v[i]);// libera as linhas da matriz 
        free (v);                       // libera a matriz 
        return (NULL);                  // retorna um ponteiro nulo 
}

float *Alocar_vetor_real (int n)
{
        float   *v; // ponteiro para o vetor

        if (n < 1) // verifica parametros recebidos
        {  
             printf ("** Erro: Parametro invalido **\n");
             return (NULL);
        }

        // aloca o vetor
        v = (float *) calloc (n+1, sizeof(float));
        if (v == NULL) 
        {
             printf ("** Erro: Memoria Insuficiente **");
             return (NULL);
        }
        return (v);    /* retorna o ponteiro para o vetor */
}

int *Alocar_vetor_int (int n)
{
        int     *v; // ponteiro para o vetor

        if (n < 1) // verifica parametros recebidos
        {  
             printf ("** Erro: Parametro invalido **\n");
             return (NULL);
        }

        // aloca o vetor
        v = (int *) calloc (n+1, sizeof(int));
        if (v == NULL) 
        {
             printf ("** Erro: Memoria Insuficiente **");
             return (NULL);
        }
        return (v);    /* retorna o ponteiro para o vetor */
}

float *Liberar_vetor_real (int n, float *v)
{
        if (v == NULL) return (NULL);
        if (n < 1) // verifica parametros recebidos
        { 
                printf ("** Erro: Parametro invalido **\n");
                return (NULL);
        }
        free(v);        // libera o vetor
        return (NULL);  // retorna o ponteiro
}
float Numero_maximo(int n,float *v)
{
	int i;
	float saida;
	saida=v[0];
	for(i=1;i<n;i++)
		if(saida<v[i])
			saida=v[i];
	return saida;			
}

